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Abstract 

One of the most recently developed ground source heat pump (GSHP) models which accounts for real-time, 
comprehensive surface energy balance phenomena, was modified in this paper to obtain the solution for a non- 
homogeneous domain. Alternating direction implicit and fully explicit finite difference methods were utilized and 
numerical results were compared to field measurements. The developed model was then used to evaluate effective- 
ness of an intermediate insulation layer obtained from tire recycling, Tire Derived Aggregate (TDA). Employment of 
TDA is demonstrated to be beneficial for increasing energy extraction rates from the ground to maximum values of 
about 17 % for cold climates, and about 5 to 6 % for moderate climate conditions. Existence of an optimum configu- 
ration for the TDA blanket, which maximizes the energy gain, also was detected. Potential energy saving for GSHP 
operation on an annual basis, with the introduction of the TDA blanket was demonstrated as practical. Validation of 
these findings with field data of longer duration are necessary and currently being planned. 

Keywords: Ground source heat pump, tire derived aggregate, alternating direction implicit, non-homogeneous, 
insulation, energy efficiency 



1. Introduction 

Ground pipes are the key unit of ground source heat pumps (GSHPs) making the link between the ground, as 
source/sink of energy, and the heat exchangers inside the heat pump [1]. Due to lack of high computational power, 
early numerical models for ground pipes were centered around line source heat theory, representing ground pipes as a 
line source in an infinite or semi-infinite medium [2]. Inherent limitations associated with line source theory, include 
the need for calculating the source strength from actual field data and its inability to include effects of circulating 
fluid. The need for overcoming these limitations led to development of a modification of the line source model 
proposed by Metz [3], where the source strength was based on the house load and the surface temperature effects were 
taken into account by introducing the Kusuda [4] temperature model. Metz also incorporated the thermal interference 
effects between nearby coils using the principle of superposition. However, the bigger challenge of assessing the 
effects of fluid characteristics on the amount of heat extraction from the ground has not yet been addressed. The next 
major advancement was by Mei [5, 6]. He proposed a new model based on energy balance between the pipe and 
surrounding soil material accounting for real-time effects of circulating fluid temperature and properties. One of the 
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main advantages of Mei's model was its ability to consider actual heat source strength by introducing time-dependent 
fluid temperatures instead of a constant source strength as in a line source model, which is not an actual representation 
of reality. The model is therefore capable of calculating solution domain temperatures in three dimensions, as well 
as accounts for the effects of different pipe material properties, circulating fluid flow rates and thermal properties, 
thermal backfill material and thermal interference between pipes [5, 7]. Mei could successfully predict the output 
water temperatures from the ground with the average deviations from the actual measurements being around 1-3 °C. 
Despite remarkable advancements in the modeling of ground pipes in these models, they still required adjustments 
to account for more elaborate surface energy balance equations used in Liston and Ling's work [8, 9]. Their work 
takes into account the effects of radiation, convection, and conduction heat transfer mechanisms. The link between 
these surface energy balance models and GSHP models was the key to a more comprehensive model for ground pipe 
analysis. These models were merged into surface energy balance equations by Demir [10] to simulate ground pipe's 
performance. The proposed model by Demir takes into account the surface energy balance equations, fed by the local 
meteorological data as upper boundary conditions and solves for the temperature distribution of the domain. Building 
on this past work, the aim of this paper is to evaluate the impact of an intermediate insulating layer on the performance 
of a ground pipe using an extended version of the Demir model. 

Tire Derived Aggregate (TDA) is the proposed insulation material to be used in conjunction with ground pipes. 
TDA mainly consists of chopped pieces of used tires in a variety of nominal sizes ranging from 1 to above 5 inches 
[11]. The idea of production of TDA from used tires, also referred to as tire chips, tire shreds, and tire mulch, was 
first initiated by Humphrey [12]. New applications for TDA, mostly in civil engineering, were proposed based on 
its unique physical characteristics. TDA's relatively low density compared to conventional backfill makes it a viable 
alternative fill material where a lighter fill material is desired in construction works [13]. Its low thermal conductivity 
on the other hand, provoked the thought of utilizing TDA as an alternative insulation material in the base insulation 
and some agricultural applications to modulate temperature fluctuation on the soil surface [12, 11, 14]. There are only 
a few studies in the literature that focus on the properties of TDA measuring physical properties [14], compaction 
densities for different sizes of chips [11] and also thermal properties of tire chip samples [11, 15, 16]. A list of thermal 
properties and density of tire shreds reported in these studies is presented in Table 1 . A thermal conductivity value of 
0.29 WmT l °C~ l , specific heat of 500 JKg l °C l and density of 720 Kgm~ 3 were chosen as representative properties 
of TDA material for this paper. 

In the present paper, the energy balance model developed by Demir was first adjusted to account for the solution 
of the governing heat equation with an intermediate insulation layer. As a result, the modified model is capable of 
solving the heat conduction problem through a non-homogeneous medium. Convergence behavior of the model was 
tested with two finite difference schemes, fully explicit and alternating direction implicit (ADI). The successful model 
was subsequently selected for the parametric study with the TDA layer. Thermal properties of TDA were introduced 
to the non-homogeneous model to assess its performance as an insulator blanket on top of installed pipes. Thickness 
and placement of the TDA layer was altered with depth to measure its impact on the results. Effectiveness of the 
TDA blanket for GSHP pipe performance was evaluated by comparing the numerical results with experimental data 
obtained from Demir's work for 37 days of climate condition in Turkey [10]. To account for different surface energy 
balance situations, several runs of the model for varying climatic conditions were simulated and results are discussed. 
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Table 1 : TDA thermal and physical properties in literature 



Reference 


Material 


Thermal conductivity 


Specific heat 


Density 






(Wm- lo C- 1 ) 




(KgrnT 3 ) 


Wappet (2006) 


Tire shred 


0.564 


507 


641 


Shao(1995) 


Tire chips 


0.149-0.164 


NA 


513 


Humphrey (2002) 


Tire chips 


0.29 


1.15 


720 


Moo-young (2003) 


Tire chips 


NA 


NA 


1060-1100 



2. Material and methods 

The surface energy balance equations are first introduced in detail. The second subsection describes the site infor- 
mation, GSHP specifications and also the way pipes are modeled from physical domain into the numerical domain. 

2.1. Surface energies 

The surface boundary condition takes into account the effects of energy balance due to variety of mechanisms 
responsible for surface-ambient heat interaction. The total energy balance on the ground surface (Q„ Watts) can be 
written as [9, 8]: 

Qt = Qc + Q e + Qh + Qie + Qu + Qst + Qp (1) 



Q c - Conduction heat flux through snow layer or ground surface (Watts) 

Q e - Turbulent exchange of latent heat (Watts) 

Qh= Turbulent exchange of sensible heat (Watts) 

Qi e = Emitted long wave radiation heat flux (Watts) 

Qu= Incoming long wave radiation (Watts) 

Q si - Solar radiation reaching the surface of earth (Watts) 

Q p - Heat flux due to precipitation (Watts) 

Heat conduction through the snow and ground layers can be written as [8]: 

Q c = -(Ts0-T b )(^r + ^r i (2) 

A 5 Kg 

where T s0 is the snow surface or ground surface temperature depending on whether the surface is covered with snow or 
not, Th is the ground temperature at the bottom of topsoil layer (thickness of first discretization element in y direction, 
Ay), z s and z g are the thicknesses of snow and the top layer of soil (Ay), respectively. A', and {WmT Xo K~ l ) are the 
thermal conductivities of the snow and soil layers respectively. 
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Turbulent exchange of sensible and latent heat, Qi, and Q e are given as [8, 9]: 

Qh=PaC p , a D h {(T a -T s0 ) (3) 

Q e = p a L s D e ^0.622 e -^^) (4) 

The exchange coefficients for sensible and latent heat, £)/, and D e , and the stability function £ are defined as 
follows: 

K 2 U 

D e = D h =-p±- 2 (5) 
(lnz/zo) 2 

^rrk (6) 

where p a is the air density equal to 1.275 kgm~ 3 , C p , a is the specific heat of air assumed to be 1000 Jkg~ l °K~ l , T a 
is the air temperature (°K), e s o is the surface vapor pressure (Pa), e a is the atmospheric vapor pressure (Pa), P a is the 
atmospheric pressure (Pa), L s is the latent heat of sublimation of snow (Jkg" 1 ), k is the VonKarman's constant, U s is 
the wind speed (ms _1 ) at reference hight z (m), and z 0 is the roughness length set to 0.5 m. 
The Richardson number, R it is defined as: 

p gz(T a - T sQ ) 

Ri = T a Ul (?) 

where g is the gravitational acceleration (ms~ 2 ), and U z is the wind speed at the elevation z relative to reference hight 
Zo- The emitted longwave radiation, Qi e , was given by: 

Qie = eMT sQ ) 4 (8) 

where e s is the soil surface emissivity assumed to be equal to 0.98 , and cr is the S te fanBoltzman constant equal to 
5.6704 • 10~ 8 (Wm- 2 K- 4 ). 

The incoming longwave radiation was given by the empirical description [17]. 

Qu = 1-08(1 -e- mu ^)o-(T a f (9) 



2353 

log w e a = 11.40-— (10) 

1 dp 

where T dp is the daily dew-point temperature (° K). The incident solar radiation reaching earth's surface can be 
described as follows: 

Q si = (1 - Albedo)[S m + S a Re(exp(iwt + 00)] (11) 

where S m is the mean annual solar radiation (Wm~ 2 ), S a is the amplitude of surface solar radiation (Wm~ 2 ), w is the 
angular velocity (rad), and (f>\ is the phase angle (rad). Heat flux due to precipitation can be written as [10]: 
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Q P = I R C p , w (T a - T sQ ) (12) 

where I R is the rain intensity (kgm~ 2 s), and C PtW is the specific heat of water assumed to be equal to 4186 Jkg~ l K~ l . 

As some of the energy equations are a function of the surface temperature of the solution domain that need to 
be numerically solved, utilization of an iterative method was necessary. The surface energy equations, which are 
obtained from the daily meteorological data, were solved for the unknown surface temperature values, T s q, using the 
Newton-Raphson method. This iterative scheme solves for temperature values on the surface using equation (13): 



f(T s0 ) 
'""/'('/„» 

where f(T s o) takes the following form: 



f(T sQ ) = Q c (T sQ ) + Q e (T sQ ) + Q h (T sQ ) + Q le (T sQ ) 

+Qu(Tso) + Q si (T s0 ) + Q p (T s0 ) = 0 (14) 

The iteration continues until the temperature difference between two consecutive iterations becomes less than 0.1 

° C. 

2.2. Ground pipe modeling 

The experimental set-up for the ground pipe (Polypropylene Random Copolymer-PPRC material) heat exchangers 
consists of three parallel horizontal pipes of total length of 40 m (L) and \ inch diameter, 3 meters apart, buried at 
the depth of 1 .8 m (j P i pe ) installed at Yildiz Technical University, Turkey. A GSHP having 4 kW heating and 2.7 
kW cooling capacity was used. Taking into account the symmetry of the problem and assuming there is no thermal 
interaction between pipes, the solution domain (Fig. 1) was considered to be from the center line of the middle pipe 
to the mid-span distance between pipes (1.5 m), in the x-direction, and from ground surface to the farfield (5 m) in 
the y-direction. Water was circulated through loops with the flow rate equal to 0.427 (m 3 h~ l ) starting on Dec 13 th 
2005. Local meteorological data were employed to force the upper boundary fluxes based on the heat flux equations 
described previously [10]. Water temperature in the inlet and outlet of the ground pipe, as well as temperature profiles 
at selected distances from the pipe were monitored for the duration of the experiment. Soil thermal conductivity 
(K s ) and thermal diffusivity (a s ) were equal to 2.18 WmT l °K~ x and 6.8 ■ 10~ 7 m 2 s~ l , respectively. To account for 
the three-dimensional behavior of the pipe and the surrounding soil, the effect of water flow was considered along 
the pipe direction. The third dimension of the problem was modeled by splitting the physical domain in the pipe 
direction to a series of cross sections (slices) of soil profiles for each time step, including the nodal temperature of 
water at the pipe's location. Fig. 2 shows how the slices are spaced in the pipe direction to cover the temperature 
distribution of the entire 3D domain. At each time step, nodal temperatures of each cross section were obtained and 
subsequently updated for the next slice along pipe's length via equation (16) to achieve a temperature distribution of 
soil and water at the end of the pipe (L). The same process was repeated for future time steps until the end of the 
simulation time. A schematic representation of the cross sections is depicted in Fig. 1. The model was run for two 
separate scenarios: a homogeneous (soil) and a non-homogeneous (with TDA blanket) medium. Each scenario used 
the same pipe configuration to evaluate effects of the TDA layer on the output water temperature, as well as the soil 
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j (TDA bottom) 
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Soil 

TDA 



Soil 



Figure 1 : Configuration of pipe, TDA and soil layers and solution domain discretization 




Figure 2: Schematic of slices of 3D domain in pipe directions 



profile. Description of numerical schemes and results are described in the following section. 
3. Numerical algorithm 

The three-dimensional temperature distribution in the soil media was modeled by solving the governing heat con- 
duction equation for the profile and incorporating the heat flow rate forced by circulating water. The temperature gra- 
dient in the pipe material is small enough to be neglected allowing the heat equation to be solved in two-dimensional 
geometry by inclusion of the water temperature in the domain [18]. The governing equation (15) was solved for the 



6 



entire domain in two-dimensions for the pure heat conduction first, and the solution was then updated along pipe's 
length to obtain the solution in three dimensions. 



1 dT(x,y,t) = d 2 T(x,y,t) | d 2 T(x,y,t) 
a dt dx 2 dy 2 

The output water temperature equation (16) along the pipe direction, calculated based on the analytical solution 
for energy balance between surrounding soil medium and pipe [19], constructs the link between the fluid and soil 
temperature in the model. 

Tf, ou , = T s -{T s -T u )e^tf (16) 

where a is the thermal diffusivity of the medium through which heat travels, Tf jOUt is the water output temperature 
coming out of a pipe of length L (m). T s is the surrounding soil temperature, Tjj is the initial water temperature 
entering pipe, and K s is the soil thermal conductivity (Wm~ lo C~ 1 ). m is the mass flow rate (Kgs~ l ) and C p j is the 
specific heat of water (Jkg~ l °C _1 ). 

Boundary and initial conditions for the solution domain are described as follows. Initial temperature distribution 
of the soil profile (T t ) was obtained from the Kusuda model, equation (17)[4]: 

Tt = T(x,y),t = 0 

^ = 0 x — x 

%=0,x = 0 

Q,(W/m 2 ),y = 0 

T(y, t) = T(y, f) Kumda , y = y max 



T(y,t) K usuda = T avg + T amp e y cos(2tt- ~yJ-^) (17) 

where T avg is the average annual surface temperature, T amp is the annual amplitude of fluctuation of surface air 
temperature, y is the depth from ground surface (m), a s is soil thermal diffusivity (m 2 s~ l ), P is the duration of a year 
in seconds, and t is the time of the year in seconds. 

Partial derivation of temperature in x-direction was written with a central difference scheme: 

ffr T n , .-2T n . + T" . 

— +0(Ax) £ 



dx 2 ,(iJ) (Ax) 2 

T" ,-2T n . + T" . x 2 T n 

(Ax) 2 " ~ (Ax) 2 

where: 



(18) 



= Ti-ij - 2T t j + T i+ ij (19) 
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Figure 3: Two dimensional grid of solution domain 



Using the same discretization for y-direction, the governing heat equation therefore can be expressed in terms of 
nodal temperature values in its explicit form as in equation (20). Next, a scheme was selected to march through time. 
As depicted in Fig. 3, the solution domain was discretized in n x + 1 nodes in x and n y + 1 nodes in the y-direction, of 
which inner domain was used to solve for the temperature distribution of each cross section of the physical domain. 
The nodes on the boundary were separated to force boundary conditions in x and y direction. 



6 2jn S 2 y T" 



(20) 



At ~ L (Ax) 2 (Ay) 2 

Two numerical solution schemes were employed to solve for temperature distribution of the solution domain. Nu- 
merical models were constructed in MATLAB. Alternating Direction Implicit (ADI) and also a fully explicit method 
were employed and the results were compared. Solution of the heat conduction equation with ADI method requires 
writing the difference formulas separately for rows and columns in two consequent time steps to solve for nodal tem- 
peratures [20]. A scheme of the two-step process for rows and columns is depicted in Fig. 4. The difference form of 
the governing equation for first and second time step then took the form of equations (21) and (22): 



> j 



T" — w"r"+l 1T»tl i f n+1 \ , r (T n OT n _i_ T n \ 

1 i,j ~ r ^ Z - L i,j + 1 i+l,j> + 'V LL i,j + 1 i,j+l> 

rrn+2 rpn+l _ t .( r rn+\ ^t^h+I , rrn+l \ , /rpn+2 ^^rn+2 , \ 

1 i,j 1 ij ~ r ^ i-lj L1 i,j + 1 i+\,j> + 'V i,7-l zi ij + 1 i,j+V 



where: 



a At 



a At 



(Ax) 2 ' ry ~ (Ay) 2 >'*-'y-' 

Rearranging these formulas to have unknowns and knowns of each time step on left and right hand side, respec- 
tively, gives: 

-rTtfj + (1 + 2r)T?f - rT^ . = rT"^ + (1 - 2r)7?. + rT n i j+x (24) 
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(21) 



(22) 



(23) 



-^£h + (1 + 2r)7** - rT™ = rT^j + d - 2r)T^ +1 + . 



r»+2 „Tti+2 _ „Tn+\ 



n+1 | „7in+l 



(25) 
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Figure 4: ADI sweeping through rows and columns 

Fully explicit finite difference formulation of heat conduction was used as second method of marching in time, 
equation (26). 



r;;; 1 - r- = r(r^. - 27- + t? +1j ) + Kr^, - 27* + r^ +1 ) 

Rearranged form of explicit equation will take the form of equation (27): 

7- +1 = rf* + rTf +1> . + rTl hl + r^ +1 + (1 - 4r)7£ 



(26) 



(27) 



Once the homogeneous case formulation was performed, the next step was adjustments to the difference equations 
to take into account the effects of the internal boundary conditions on the top and bottom of the TDA layer. To meet 
the temperature and flux conditions on internal boundaries, the energy balance on the interfacial nodes of TDA and 
soil was obtained by writing the nodal fluxes and the energy storage term for the control volume around each node. 
Fig. 5 depicts the energy balance for a node on the top and bottom of TDA layer interface, the control volume around 
the node, and the heat fluxes involved. The resulting energy balance is shown in equation (28): 



qi + q 3 - q 2 - qn = (pC p ) avg ■ Ax ■ Ay 



dt ' 



where: 



/ \ _ (PC,,)MPC,,)TD,\ 

\P^p)avg — 2 



K, 



avg 



and the values of heat fluxes are: 



(28) 



i-ij 



Ax 



■Ay 



(29) 



qi - K av 



i+lj 



•J 



Ax 



■Ay 



(30) 
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Figure 5: TDA layer top and bottom interfaces energy balance 



T7i y« 

93 = ■ - ' • Av (31) 

Ay 

Y*n r p n 

94 = K TDA ■ — ■ Ax (32) 

Ay 

Substituting the flux equations into the energy balance equation and rearranging the equations results in the intro- 
duction of new parameters in the difference equation as described below. 



_ K avg At _ K TDA At = K s At 

ravg ~ (pC p ) avg ' (Ax) 2 ' VtDA ~ (pC p ) avg ' (Ax) 2 ' r ° ~ (pC p ) avg ' (Ax) 2 

The difference formulas obtained previously for the ADI method were revisited, resulting in new difference formu- 
las for top and bottom interfacial nodes of TDA and soil layer. New interfacial formulas for rows, equations (34) and 
(35), as well as for columns, equations (36) and (37) were obtained. Similar changes applied to the explicit method 
resulted in equations (38) and (39) to account for the heat flux through the TDA top and bottom layers. 

For TDA top row: 

-ranTt-tj + (1 + 2r avg )T^ 1 - r avg T^ = r s T^ + (1 - r, - r TDA )T*j + r T T; . +l (34) 
for TDA bottom row: 

-r mg T?lj + (1 + 2r avg )T?f - r^T^ . = r T T^_, + (1 - r, - r TDA )T»j + r s T? j+l (35) 
for TDA top column: 

-r s T?fi + (1 + r T + r s )T% 2 - r T T^ x = r avg T^j + (1 - 2r avg )T?f + r avg T^j (36) 
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for TDA bottom column: 



-r T T^ + (l + r T+ r s )T?f - r s T^ = r mg T^j + (1 - 2r flvg )7- +1 + r avg T^. (37) 
for explicit TDA top layer: 

T" +1 = r avg T»_ Uj + r avg T» +1J + r s T» M + r T T» j+1 + (1 - 2r flvg -r s - r r )T£ (38) 
for explicit TDA bottom layer: 

- r av Jl x . + r mg T? +1J + r T T»._ x + r s T» j+l+ (l- 2r avg -r s - r T )T^ (39) 

4. Results and discussion 

Based on Demir's field data, the ADI method first was used successfully to model the exact same problem solved 
for the homogeneous case with soil material [10]. The homogeneous case was then solved with fully explicit methods 
to compare the results with ADI. The results of the ADI method output water temperature are plotted in Fig. 6 versus 
measured input/output water temperature to/from ground during the thirty seven days of the field experiment. Two 
mesh sizes, 0.05 and 0.1 m, were used in ADI runs to check the sensitivity of the results to space discretization, while 
time increments were held constant at 1800 seconds for all runs. The modeled output temperatures for both mesh sizes 
were in good agreement with the measured data for almost all the daily data values. Similarly, a plot of numerical 
results were obtained for the explicit model as depicted in Fig. 7. Mesh size 0.1m and time increment of 1800s were 
used for explicit runs. The explicit method results showed fairly better agreement with the field measurements. 



i i 

Input water Temp. 

— a — Output Water Temp(Experiment) 
— e— Output Water Temp(ADI) dx=0.1 
— * — Output Water Temp(ADI) dx=0.rj5 
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Figure 6: Water temperatures to/from ground, experimental data vs. ADI method results 
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Figure 7: Water temperatures to/from ground, experimental data vs. explicit method results 



After gaining the confidence in the model's preciseness for homogeneous runs, the two methods were examined 
for non-homogeneous case by the introduction of TDA blanket as an internal layer above the pipe. Mathematical 
derivation of the stability analysis, developed by Gao [20], shows that the ADI method is unconditionally stable for 
homogeneous, two-dimensional heat diffusion problems . It is well established in the literature that the ADI scheme 
is fairly efficient as it requires the solution to a tri-diagonal matrix to obtain nodal values. However, as demonstrated 
by Rivera [21], the Von Neumann stability analysis method used for numerical boundary conditions yields only a 
necessary condition for stability because it does not consider the overall effect of the boundary conditions between 
subdomains. Rivera explains the stability problem associated with the solution of PDEs when the original domain is 
decomposed into overlapping subdomains. He also states that explicit predictor method can be used for interfacial 
nodes to obtain stable solutions. The same stability analysis was used by Gao [20] to prove the stability of the 
ADI method, therefore it seemed reasonable to be more cautious about the non-homogeneous case results with the 
ADI method. However, a successful example of ADI method used for a non-homogeneous case has reported good 
performance except over small times and large distances wherein the numerical round-off error causes considerable 
deviation from the analytical solution [22]. Moreover, Trefethen [23] showed that the stability of finite difference 
models of hyperbolic initial-boundary value problems can be highly affected by boundaries when there are additional 
boundaries or internal interfaces in the physical model. 

Considering the history of the stability conditions for the case of internal boundary conditions, the codes were 
modified to account for the TDA blanket as prescribed by the energy balance equation (28). The ADI method was 
found to be unstable for the non-homogeneous case with Ax-Ay-Q. 1 m, and At- 1800 s. Explicit method results with 
the same time and space discretization, on the other hand, were stable and successfully used for a parametric study for 
the problem with the TDA blanket. The explicit method stability criteria, 1 - 4r > 0 and 1 - 2r avg - r$ - fr > 0, were 
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met for both homogeneous and non-homogeneous cases. Several time and space discretizations were used to assess 
sensitivity of the ADI method but none resulted in a stable solution of the problem. In a diagnostic attempt, thermal 
properties of the intermediate layer were gradually changed with respect to the soil properties to explore the reasons 
behind the instability of the ADI method. Starting from the homogeneous case, the specific heat of the intermediate 
layer was lowered to 5, 20, 30, 50 and 80 percent of the soil specific heat value to assess the sensitivity of the results 
to changes in intermediate layer's properties. The output temperature values exhibited dramatic increase compared 
to the measured data for the homogeneous case for more than 5 percent change in the properties of the intermediate 
layer. One reason for such behavior of the model lies behind the new parameters emerged from the energy balance 
equation at the interfaces; r avg , rjDA, r s . The values of these parameters, calculated for TDA properties, are one order 
of magnitude different than the values for the same nodes in the homogeneous case. This considerable difference 
probably caused the sweeping process through the columns more erroneous as the nodal values at the interface are 
significantly different from the adjacent top and bottom nodes. This assumption was verified by solving for the entire 
domain with a semi-implicit method sweeping only in row direction where consistently similar results were obtained 
compared to the fully explicit method. The interfacial parameters for ADI may also introduce round-off errors that get 
magnified in the iterative solution of the upper boundary condition and subsequently exhibit nonconvergent behavior 
over time. 

A parametric study was designed to evaluate the effects of thickness, and configuration of TDA on simulation 
results, as well as its dependence on different meteorological conditions. In the first round of runs for the explicit 
model with TDA, the top position of TDA was fixed at 0.5 m and three thicknesses of TDA equal to 0.2, 0.5 and 1 
m were employed. Selected temperature snapshots of solution domain were plotted in selected horizontal distances 
from the pipe, 0.2, 0.7 and 1 .2 m to give an idea of how effective the blanket can be in winter time. Results for the 
three thicknesses of TDA are presented in Figs. 8 , 9 and 10. It was evident that the TDA blanket kept the ground 
warmer compared to the homogeneous case with the soil for almost the entire depth of the domain. The effect of pipe's 
temperature on the surrounding soil temperature was also obvious from the figures where the profile inclined towards 
the water temperature's vicinity. In all three cases, the temperature difference between the case with and without TDA 
was lower at a closer distance to the pipe than the rest of profile, even though the magnitudes varied. The interesting 
observation from the profile snapshots was that with increase in the TDA thickness while fixing pipe's location (1.8 
m) and TDA top position, the soil profile did not necessarily remain warmer with more TDA material in the ground. 
In other words, if we were to calculate the area between temperature profiles with and without TDA, there would exist 
an optimal configuration of TDA that would keep the ground warmer with TDA than with only soil. Analyzing the 
profiles showed the potential advantage of the TDA blanket to increase the flux to the pipe and subsequently increase 
the temperature of the water coming out of the ground, so the next step was quantifying this potential benefit. 

A new set of runs of the model were designed to assess the effect of different configurations of the TDA blanket 
and air temperature fluctuations. It was assumed that a fixed 0.5 m thick TDA blanket was placed at 0.1, 0.3, 0.5, 
0.7, 1, and 1.2 m (top position) from soil surface to examine the effects on temperature of water coming out of the 
ground. A single aggregate temperature value, equal to the average output water temperature over thirty seven days of 
measurements (r™ A and T*°' g l ) was assumed to be representative of the effectiveness of TDA versus soil for the entire 
period of simulation. The difference between the averaged temperatures [A(r™ A - T^™ 1 )] was then used to calculate 
the increase in earth energy absorption rate provided by employment of the TDA blanket. Increases in the energy 
absorption rates versus the location of the TDA's top position are plotted in Fig. 11. This procedure was repeated for 
three separate mean annual surface temperatures (T avg ) and fluctuation amplitudes (T amp ), representing mild (lower 
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Figure 8: Soil temperature profile at selected distances from pipe's centerline for 0.2 m thick TDA blanket 
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Figure 9: Soil temperature profile at selected distances from pipe's centerline for 0.5 m thick TDA blanket 



fluctuation amplitudes, 10.8 °C) to cold (higher fluctuation amplitudes, 15 and 25 °C) weather conditions, to somehow 
quantify the effects of TDA blanket as compared to the soil, in different hypothetical climates. Results, depicted in 
Fig. 11, showed that increasing energy absorption rates with increasing depth of the TDA blanket reaches a maximum 
level. Placing the TDA below that depth decreases the absorption rate. The other trend observed was that the higher 
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Figure 10: Soil temperature profile at selected distances from pipe's centerline for 1 m thick TDA blanket 
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Figure 1 1 : Increase in average absorbed energy from the ground with TDA blanket for different weather scenarios 



the fluctuation of the surface air temperature, the higher the effectiveness of the TDA in terms of increasing the energy 
absorption compared to the soil profile. Increases in energy absorption rates for high temperature amplitudes (15 and 
25 °C) are considerably higher than the mild temperature fluctuation amplitude (10.8 °C), where mild climate is the 
natural condition of the site. The corresponding percentage increase in absorbed energy from the ground with varying 
burial depth of the TDA is presented in Fig. 12. Results suggest that an increase in absorbed energy can be as high as 
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Figure 12: Percent increase in absorbed energy from the ground 
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Figure 13: Maximum output water temperature difference (Degrees C) for different weather scenarios in 37 days of experiment 



around seventeen percent for relatively cold climates to as low as five to six percent for fairly mild climates. 

A similar analogy was drawn by Incropera where he discussed the possible existence of an optimum insulation 
thickness for radial systems involving the flow of a refrigerant and an insulation medium around the conducting 
tube [24]. The problem definition was centered around the idea of competing effects associated with an increase 
in the thickness of insulator medium, where heat conduction and convection both existed. The problem was then 
solved for steady-state. It was acknowledged that "although the conduction resistance increases with the addition of 
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insulation, the convection decreases due to increasing outer surface area". Followed by selection of a series of different 
thicknesses for the insulator, it was revealed that there exists a "critical radius" that maximizes total resistance to heat 
transfer. Even though this was not the perfect match for our problem with respect to boundary conditions and time- 
dependency, Incropera's example suggested the existence of an optimal placement for insulators which was observed 
in our problem as well. 

It was noteworthy to check the maximum temperature differences attainable with the introduction of TDA. Maxi- 
mum temperature data are important as several GSHP designers use extreme annual temperature values to determine 
pipe length and pump sizing. There are new approaches in the literature compared to peak load design criteria which 
take into account the real time performance of a GSHP system [25, 26]. Assuming that design guidelines will be 
based more on real-time variation of air and soil temperature in the near future, an emphasis that has already started 
in research studies [27, 28, 26, 29], it was essential to monitor maximum output water temperature differences for this 
problem. The maximum temperature differences for three climatic conditions are plotted in Fig. 13. TDA demon- 
strated to be able to increase output water temperature from about 0.4 °C to as high as 1 °C, an increase of temperature 
that might have a meaningful impact on the performance of the heat pump if it occurs frequently throughout the year. 
A comprehensive analysis of the heat pump's performance for a longer time scale including all other units of a GSHP 
(e.g. heat exchangers, compressors, and condensers), might reveal advantages of TDA more clearly in terms of long 
term energy savings. Research works have emphasized that monitoring and modeling performance of GSHPs for 
small time steps including interaction of internal units of heat pumps are inevitable for component optimization of 
these systems [29, 30, 31] 

5. Conclusions 

Performance of a GSHP's pipe was investigated for thirty seven days of winter with the introduction of a new 
insulation material, Tire Derived Aggregate (TDA). Demir's experimental data were used and the code was modified 
to account for the intermediate layer [10]. Two numerical solution schemes were successfully employed to solve for 
the homogeneous case, whereas only the fully explicit method yielded stable results for the non-homogeneous case 
with the TDA blanket. Here is a summary of key results and findings: 

• ADI and explicit methods worked equally well to model homogeneous heat conduction through soil 

• The fully explicit method provided reasonable, stable results with the introduction of the intermediate TDA 
blanket as an insulator 

• The TDA layer proved effective in increasing the amount of energy absorbed from earth by as much as seventeen 
percent 

• On an average temperature basis, the TDA blanket proved to have a relatively higher performance in colder 
climates than mild ones 

• The existence of an optimum configuration for TDA blanket was detected 

• Potential energy saving of GSHP performance with the TDA blanket needs to be investigated with more field 
data and for a longer periods of time 
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